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ABSTRACT 

We present full-orbit phase curve observations of the eccentric (e ~ 0.08) transiting hot Jupiter 
WASP-14b obtained in the 3.6 and 4.5 /iin bands using the Spitzer Space Telescope. We use two 
different methods for removing the intrapixel sensitivity effect and compare their efficacy in decou¬ 
pling the instrumental noise. Our measured secondary eclipse depths of 0.1882% ± 0.0048% and 
0.2247% ± 0.0086% at 3.6 and 4.5 ^m, respectively, are both consistent with a blackbody tempera¬ 
ture of 2402 ± 35 K. We place a 2cr upper limit on the nightside flux at 3.6 /rm and find it to be 
9% ± 1% of the dayside flux, corresponding to a brightness temperature of 1079 K. At 4.5 /rm, the 
minimum planet flux is 30% ± 5% of the maximum flux, corresponding to a brightness temperature 
of 1380 ± 65 K. We compare our measured phase curves to the predictions of one-dimensional radia¬ 
tive transfer and three-dimensional general circulation models. We find that WASP-14b’s measured 
dayside emission is consistent with a model atmosphere with equilibrium chemistry and a moderate 
temperature inversion. These same models tend to over-predict the nightside emission at 3.6 pern, 
while under-predicting the nightside emission at 4.5 /rm. We propose that this discrepancy might be 
explained by an enhanced global C/0 ratio. In addition, we find that the phase curves of WASP-14b 
(7.8 Mjup) are consistent with a much lower albedo than those of other Jovian mass planets with 
thermal phase curve measurements, suggesting that it may be emitting detectable heat from the deep 
atmosphere or interior processes. 


1. INTRODUCTION 

Over the past two decades, observations of exoplanets 
have uncovered a stunning diversity of systems. Major 
improvements in the capabilities of space- and ground- 
based telescopes in recent years have led to the discovery 
and characterization of hundreds of new exoplanets, cov¬ 
ering a broad range of orbital properties, inter ior and at- 


mospheric compositions, and host star types (Han et al. 
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20141. Meanwhile, these same technological advances 
have enabled the detailed study of the atmospheric prop¬ 
erties of the brightest and largest planets through high- 
precision photometry and spectroscopy. Most of these 
targets belong to the class of gas giant planets known as 
hot Jupiters. The high levels of incident flux, slow ro¬ 
tation, and potentially large temperature gradients be¬ 
tween hemispheres characteristic of these planets allow us 
to test atmospheric models in a new regime unlike any 
found in the Solar System. In addition, hot Jupiters are 
some of the most favorable targets for measuring elemen¬ 
tal abundances, since most material is not in a condense d 
form at these high temperatures (e.g.. Line et ^|2014 1. 

Atmospheric circulation models of hot Jupiters pre¬ 
dict broad super-rotating equatorial jets that circulate 
energy between the day and night sides, with the pre¬ 
cise effect of these winds on the day-to-night tempera¬ 
ture contrast being strongly dependent on the particu- 
lar orbital and atmosph eric properties of the planet (see 
Heng fc Showman|2015 and references therein). As a re- 


sult of their short orbital periods, hot Jupiters have high 
atmospheric temperatures and emit relatively strongly 
at infrared wavelengths, allowing for direct measurement 
of their atmospheric brightness as a function of orbital 
phase. These phase curves can th en be converted into a 
longitudinal temperature profile (|Cowan & Agol||2008|). 
By comparing the measured phase curves to theoreti- 
cal phase curves generated by atmospheric models, we 
can constrain fundamental properties of the atmosphere, 
such as the efficiency of heat transport from the day side 
to the night side, radiative time scales, wind speeds, and 
compositional gradients between the day and night sides. 

To date, well-characterized phase curve observations 
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have been published for eleven planets: v And b (Har- 


jKnutson et al.| 

son et al.||2UU9 


Mr? 


3) 


rington et al. 2006 Crossfield et al. 

2009al|20l4 

Hb 80606b jLaughlin et al.| 


20101, HD 18973M 


HUT49026b (Knut- 


HAi’-H-Tb (Horucki et al. 


Kepler-7b (Demory et al. 


Welsh et al. 
CoRoi'-ib (ISnellenl 


mi: 


MTU), 


120131), HA\''-P-2b ( 

Lewis et al. 120131, and 

Crosslield et al. 

20 

12 Zellem et al. 

2014). 


The majority ot these observations were carried out using 
the Spitzer Space Telescope while the rest were obtained 
at optical wavelengths by the CoRoT and Kepler mis¬ 
sions. Recently, the first spectroscopic phase curve was 
obtained for the hot Jupiter WASP-43b usin g the Hub¬ 


ble Space T elescope between 1.2 and 1.6 ^m (Stevenson 


Die space I t 

eTidllMI ). 

in this pai 


paper, we present full-orbit phase curves of the 
hot Jupiter WASP-14b in the 3.6 and 4.5 /im bands 
obtained with the Spitzer Space Telescope. Photomet¬ 
ric and radial-velocity observations of WASP-14b indi¬ 
cate a mass of Mp = 7.3 ± 0.5 Mjup and a radius of 
Rp = 1.28 ± 0 .08 Rjup, corresp onding to a density of 


p = 4.6 g/cm (Joshi et al.||2009). The planet lies on an 
eccentric orbit (e = 0.0822 ± 0.003; Knutson et al.||2014 1 
around a young host star (age ~ 0.5 — 1.0 Gyr, spectral 
type F5, M* = 1.21 ± 0.13 Mq, i?* = 1.31 ± 0.07 


Tt = 6462 ±75 K and log g = 4.29 ± 0.04; Joshi et ah] 
2009||Torres et al.||2012|), with a period of 2.24 days and 
an orbital semi-major axis of a = 0.036 ± 0.01 AU. 

The equilibrium temperature of WASP-14b is rela¬ 
tively high (Teq = 1866 K, assu ming zero albed o and 
reemission from the entire surface; |Joshi et al.|2009|) , sug¬ 
gesting that the thermal evolution ot WASP-14b rnay be 
signihcantly affected by Ohmic dissipation in a partially- 


ionized atmosphere (e.g., Batygin et al. 2011 2013 |Perna 


et al. 112010 but for debate seelRogers h Showman||2014[ 
Rogers & komacek||2014|). |Blecic et al.| p013|) analyzed 
secondary eclipse observations in the 3.6, 4.5, and 8.0/rm 
Spitzer bands and found no evidence for a thermal inver¬ 
sion in the dayside atmosphere, while concluding that 
the data are consistent with relatively poor energy re¬ 
distrib ution from dayside to nightside. ICowan fc Agol| 
(j2011|) find that, based on the dayside fluxes, the most 
highly-irradiated planets have s ystematically less effi¬ 


cient day/n ight heat circulation. Perez-Becker & Show- 
( 2013[) reached a similar conclusion from a compari- 


man 


son o: 


the fractional day-night flux differences of planets 
for which phase curves have been obtained. This over- 
all trend has been explained by hydrodynamical m odels 
( Perna et al.||2012 Perez-Becker fc Showman||2013 ). 

Phase curve observations at more than one wave¬ 
length provide complementary information about the 
properties of the planet’s atmosphere, as different wave¬ 
lengths probe different pressure levels within the atmo¬ 
sphere. Multiband measurements of the planet’s bright¬ 
ness can also be transformed into low-resolution day- 
side and nightside emission spectra, which can reveal 
differences in atmospheric composition. Only five sys¬ 
tems have phase curve o bservations at more than one 
wa velength: HD 18973 3 (Knutson et al. 2012), WASP - 
12 dCowan et al.||2012|), WASP -18 (iMaxted et a l.||2013|), 
HAT-P-2 ([Lew is et al.||2013|), and WAyP-43 (pteven- 

■’iAi /it — „i'„ __ 


son et al. 2014[). W hile single-wavelength phase curves 


can be reasonably well-matched by standard atmospheric 
circulation models, none of the multi-wavelength phase 
curves are satisfactorily reproduced by these same mod¬ 
els. This suggests that our understanding of the physical 
and chemical processes that drive atmospheric circula¬ 
tion is still incomplete. 

The paper is organized as follows: The observations 
and data reduction methodology are described in Sec¬ 
tion In Section we discuss the phase curve model 
used in our analysis and present the best-ht parameters. 
We then use our results to obtain updated orbital pa¬ 
rameters and discuss the implications of our phase curve 
fits for the planet’s atmospheric dynamics in Section]^ 

2. SPITZER OBSERVATIONS AND PHOTOMETRY 

We observed two full orbits of WASP-14b in the 3.6 and 
4.5 Atm channels of the Infrared Array Camera (IRAC; 
Fazio et al. 2004) on the Spitzer Space Telescope. The 
observation periods were UT 2012 April 15—17 and UT 
2012 April 24—26 for the 3.6 and 4.5 /im bandpasses, 
respectively. Both observations lasted approximately 64 
hours and were carried out in subarray mode, which gen¬ 
erated 32 X 32 pixel (39"x39") images with 2.0 s inte¬ 
gration times, resulting in a total of 113,408 images in 
each bandpass. Due to long-term drift of the telescope 
pointing, the telescope was repositioned approximately 
every 12 hours in order to re-center the target, leading 
to four breaks in each phase curve observation with a 
combined duration of about 16 minutes. The telescope 
repositioning produced offsets of up to 0.2 pixels in the 
star’s position on the array after each break (Figures 
andj^. 

We extract photometry using methods described in 
previous analyses of post-cryogenic Sp itzer data (e.g., 
Todorov et al. 2013 Wong et al. 2014). The raw data 
hies are hrst dark-subtracted, nat-helded, linearized, 
and flux-calibrated using version S19.1.0 of the IRAC 
pipeline. The exported data comprise a set of 1772 
FITS files, each containing 64 images and a UTC-based 
Barycentric Julian Date (BJDutc) time stamp desig¬ 
nating the start of the first image. For each image, 
we calculate the BJDutc at mid-exposure by assuming 
uniform spacing and using the start and end times of 
each 64-image series as defined by the AINTBEG and 
ATIMEEND header keywords. 

When estimating the sky background in each image, we 
avoid contamination from the wings of the star’s point- 
spread function (PSF) by excluding pixels within a ra¬ 
dius of 15 pixels from the center of the image, as well 
as the 13th-16th rows and the 14th and 15th columns, 
where the stellar PSF extends close to the edge of the ar¬ 
ray. In addition, we exclude the top (32nd) row of pixels, 
which have values that are consistently lower than those 
from other pixels in the array. We take the remaining set 
of pixels and iteratively trim values that lie more than 
3(7 from the median. We then calculate the average sky 
background across the image by fitting a Gaussian func¬ 
tion to the histogram of the remaining pixel values. After 
subtracting the sky background, any remaining transient 
“hot pixels” in each set of 64 images varying by more 
than 3(7 from the median pixel value are replaced by the 
median pixel value. In both bandpasses, the average per¬ 
centage of replaced pixels is less than 0.35%. 

To determine the position of the star on the array in 
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each image, we calculate the 
a circular region of radius tq 
timated po sition of the star ( 
et al. 2008 1 . We then estimate 
by computing the noise pixel 
which is defined in Section 2.2 
strument handbook as 


flux-weighted centroid for 
pixels centered o n the es- 
see, for example, Knutson 
the width of the star’s PSh' 


parameter (Mighell 2005), 
.2 of the Spuzer/JUJnjm- 


_ (E^ 


( 1 ) 


where li is the intensity detected in the ith pixel. We 
define the parameter ri to be the radius of the circular 
aperture used to calculate /3. 

We calculate the flux of the stellar target in each image 
using circular aperture photometry. We generate two sets 
of apertures: The first set uses a fixed aperture with radii 
ranging from 1.5 to 3.0 pixels in 0.1-pixel steps and from 
3.0 to 5.0 pixels in 0.5-pixel steps. The second set utilizes 
a time-varying radius that is related to the square-root of 


the noise pixel parameter jS by either a constant sc aling 
factor or a constant shift (see Lewis et al. 2013 for a 
full discussion of the noise-pixel-based aperture). The 
optimal choice of aperture photometry is determined by 
selecting data from 8000 images spanning the planetary 
transit with a total duration of 4.4 hours and calculating 
the photometric series for each choice of aperture. We 
then fit each photo metri c series with our transit light 
curve model (Section |3.1| ), compute the RMS scatter in 
the resultant residuals binned in five-minute intervals, 
and choose the values of rp and ri as well as aperture 
type that give the minimum scatter. In these fits, we 
fix the planet’s orbital parameters to the most recent 
values in the literature: P = 2.2437661 d, i = 84.32°, 


a/R* = 5.93, e = 0.0822, and to = 251.67° (Joshi et al. 
Knutson et al.|[2014). For the 3.6 fim data set, we 


2009 


hnd that a fixed aperture with a radius of 1.8 pixels and 
ro = 3.0 produce the minimum scatter. When using a 
fixed aperture, the noise pixel parameter is not needed. 
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Fig. 1.— Measured stellar x centroids (top panel), y centroids 
(upper middle panel), and noise pixel values (lower middle panel) 
as a function of orbital phase relative to transit for the 3.6 fim phase 
curve observation. The bottom panel shows the raw photometric 
series with hot pixels excised. The data are binned in 2-minute 
intervals. 
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Fig. 2.— Measured star centroids, noise pixel values, and raw 
photometric series for the 4.5 ym phase curve observations; see 
Figure^ for a complete description. 

so ri is undefined. In the 4.5 fim bandpass, we prefer a 
fixed aperture with a radius of 2.9 pixels and tq = 3.5. 

Prior to fitting the selected photometric series with 
our full light curve model, we use a moving median fil¬ 
ter to iteratively remove points with measured fluxes, x 

positions, y positions, or values that vary by more 
than 3a from the corresponding median values in the 
adjacent 64 frames in the time series. Choosing a larger 
or smaller interval for computing the median values does 
not significantly affect the number of excised points. The 
percentages of excised points are 1 . 8 % and 1 . 6 % in the 
3.6 and 4.5 ym bandpasses, respectively. 

3. DATA ANALYSIS 

3.1. Transit and secondary eclipse model 

Each full-orbit observation contains one transit and 
two secondary ec lipses. We model thes e events using 
the formalism of Mandel & Agol| (|2002 ). The transit 
light curve includes tour tree parameters: the scaled or¬ 
bital semi-major axis a/i?*, the inclination i, the cen¬ 
ter of transit time tx, and the planet-star radius ratio 
Rp/R^,, which is the square root of the relative tran¬ 
sit depth. Each secondary eclipse event is defined by a 
center of eclipse time Ie and a relative eclipse depth d, 
measured with respect to the value of the phase curve at 
mid-transit, in addition to a/i?* and i, thus yielding four 
additional free parameters: Iei, tE 2 , di, and ^ 2 - We en¬ 
sure continuity between the phase curve and secondary 
eclipse light curves by scaling the amplitudes of eclipse 
ingress and egress (when the planet is partially occulted 
by the star) appropriately to match the out-of-eclipse 
phase curve values at the start and end of the eclipse. 
The host star WASP-14 has an effective temperature 
T* = 6462 ± 75 K, a specific gravity logg = 4.29 ± 0.04, 


and a metallicity of [Fe/H] = —0.13± 0.08 (Torres et al. 
2012). We model the limb-darkening in each bandpass 


using a four-parameter non-linear limb-darkening law 
with p arameter values calculated as described in Sing| 
(2010) for a 6500 K star with log 5 = 4.50 and [Fe/HJ 
E^IO: Cl - C 4 = [-0.0192,0.7960,-0.8558,0.2983] at 
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3.6 /xm and ci - C 4 = [0.0225,0.3828,-0.2748,-0.0522] 
at 4.5 /imQ 


3.2. Phase curve model 

WASP-14b has a relatively low eccentricity of 0.08, and 
therefore the variation in the planet’s apparent bright¬ 
ness throughout an orbit can be modeled to first order 
as a simple sinus oidal function of the true anomaly / 
( Lewis et al.|2013 ), analogous in form to a simple sine or 
cosine ot the orbit al phase angle that is used in the case 
of a circular orbit ( Cowan fc Agol|[^008| ) : 

( 2 ) 


F{t) = Fo + ci cos(/(t) - C 2 ), 

where Fq is the star’s flux (assumed to be constant and 
normalized to one), Ci is the amplitude of the phase vari¬ 
ations, and C 2 represents the lag between the peak of 
the planet’s temperature and the time of maximum inci¬ 
dent stellar flux due to the finite atmospheric radiative 
timescale of the planet. Here, ci and C 2 are free pa¬ 
rameters that are computed in our fits. We also exper¬ 
imented with other functional forms of the phase curve 
that included higher harmonics, but all of them resulted 
in higher v alues of the Bayesian Information Criterion 
(see Section 3.4 for more information). 


3.3. Correction for intrapixel sensitivity variations 

Photometric data obtained using Spitzer/IRAC in the 
3.6 and 4.5 fim bandpasses exhibit a well-studied in- 
strumental effect due to in trapixel sensitivity variations 
(Charbonneau et al. 2005). Small changes in the tele- 
scope pointing during observation cause variations in the 
measured flux from the target, resulting in a character¬ 
istic sawtooth pattern in the raw extracted photometric 
series. In our analysis, we decorrelate this instrumental 
systematic in two ways. 

Our first approach to removing t he intrapixel sensitiv¬ 


ity effect is called pixel mapping (Ballard et al. 2010 
Lewis et al. 2013). In an image j, the location of the 


target on the array is given by the measured centroid 
position (xj^yj), and the sensitivity of the pixel at that 
location is determined by comparing other images with 
measured centroid positions near {xj,yj). The effective 
pixel sensitivity at a given position is calculated as fol¬ 
lows: 


Fmeas., = F, X (3) 


2 = 0 


X e 




Here, Fmeas.j is the flux measured in the jth image and 
Fj is the intrinsic flux; Xj, pj, and (3 are the measured 
X position, y position, and noise pixel parameter values. 
The quantities o '^ . are the standard 

deviations of x, y, and over the full range in i. For 
each image j, the summation in Eq. (|^ runs over the 
nearest to = 50 neighbors of the stellar target, where we 


^ Tables of limb-darkening parameter values, calculated in the 
Spitzer bandpasses for various stellar temperatures, specific grav¬ 
ities, and metallicities, can be found on David Sing’s website; 
WWW. astro. ex. ac. uk/people / sing 



Fig. 3.— Plot comparing the binned residual RMS resulting from 
the fit of the 3.6 pm photometric series to the light curve model 
using various types of instrumental noise decorrelation: pixel map¬ 
ping (solid black line) and pixel level decorrelation with no binning 
(solid red line), 16-point bins (solid green line), and 256-point bins 
(solid blue line). For comparison, the Ijy/n dependence of white 
noise on bin size is shown by the dashed black line;the trend is nor¬ 
malized at bin size n = 1 to the photon noise level corresponding 
to the median photon count over the 3.6 ^m photometric series, 
with sky background included. 


define distance as 

dlj = {Xi - Xjf + {pi - Pjf + ( 


( 4 ) 


This method in effect adaptively smoothes the raw pixel 
map, allowing for a finer spatial scale in regions where 
the density of points is high while using a coarser spatial 
scale in regions with sparser sampling. We chose this 
number of neighbors to be large enough to adequately 
map the pixel response while maintaining a reas onably 


low computational overhead (Lewis et al. 2013). Sev- 


son et al. 

2012 

Zellem et al. 

noise pixe 

parameter term in 


our WASP-14b data, including the noise pixel parameter 
term produces ^ 5 — 10% smaller residual scatter from 
our best-fit light curve solution in both bandpasses. 

The second approach is a recently proposed t echnique 


known as pixel-level decorrelation (PLD; see Deming 


et al. 2015 for a complete description). Unlike most other 


treatments of the intrapixel sensitivity effect, PLD does 
not attempt to relate the variations in the calculated po¬ 
sition of the target on the pixel to the apparent intensity 
fluctuations. Rather, it utilizes the actual measured in¬ 
tensities of the individual pixels spanning the stellar PSF 
to provide an expression of the total measured flux. We 
consider pixels lying in a 3 x 3 box centered on the star, 
which have pixel intensities Pk{t), k = 1,...9. We di¬ 
vide each 3x3 pixel box by the summed flux over all 
nine pixels in order to remove (at least to zeroth order) 
any astrophysical flux variations, giving the following re¬ 
lation: 
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The intrapixel sensitivity effect is modeled as a linear 
combination of the arrays Pk, and thus the total mea¬ 
sured intensity S is given by 


S{t) = Y. ]+btt + F{t) + h, 


( 6 ) 


\.k=l 


where F(t) is the astrophysical model, comprising the 
phase curve, transit, and eclipses. The parameters bk are 
the linear coefficients that are determined through least- 
squares fitting, and /i is a free normalization parameter 
that corrects for the overall numerical offset introduced 
by the linear s um o f pixel intensity arrays. Following 
Deming et al. (20151, we also include a linear ramp in 
time with a slope paramet er bt- 


In Deming et al. ([2015 ), PLD was applied to fitting 


Spitzer secondary eclipses of several exoplanets, includ¬ 
ing WASP-14b. PLD was found to be generally more 
effective in removing time-correlated (i.e., red) noise, re¬ 
sulting in lower residual scatter from the best-fit solution 
when compared with other decorrelation techniques. The 
best results were obtained when the photometric time 
series were binned, pixel by pixel, prior to PLD fitting, 
since binning improves the precision of the measured in¬ 
tensities in pixels at the edge of the stellar PSF and can 
be adjusted to reduce the noise on the timescale of inter¬ 
est. When fitting our full-orbit WASP-14b photometric 
series using the PLD technique, we experimented with 
fitting either unbinned or binned data, with bin sizes 
equal to powers of two up to 256 (~8.5 minutes). We 
found that larger bin sizes, comparable to or exceed¬ 
ing the occultation ingress/egress timescale of roughly 20 
minutes, cause excessive loss of temporal resolution and 
yielded unsatisfactory secondary eclipse and transit light 
curve fits. After fitting the model light curve to the data 
using PLD, we subtract the best-fit solution from the 
raw unbinned data to produce the residual series, which 
we use to evaluate the relative amount of time-correlated 
noise remaining in the data. 

In Figure!^ we compare the noise properties of each 
version of the 3.6 fim residual time series to the ideal 

\fn scaling we would expect for the case of independent 
(i.e., “white” noise) Gaussian measurement errors, where 
n is the bin size. The white noise trend is normalized at 
bin size n = 1 to the photon noise level corresponding 
to the median photon count over the observation data 
set (with sky background included). Comparing the var¬ 
ious PLD fits, we find that at small bin sizes, unbinned 
PLD gives the lowest residual scatter, while at larger 
bin sizes approaching the duration of eclipse ingress or 
egress, PLD with larger bins results in lower residual 
scatt er. The same trend is seen when fitting the 4.5 /im 
data. Deming et al. ( 2015|) found that the optimal PLD 
performance is achieved when the range of star positions 
is lower than 0.2 pixels. The range of pixel motion in 
our full-orbit phase curve observations modestly exceeds 
this limit, and we indeed find that the residual scatter 
from the pixel mapping fits is smaller than the scatter 
from any of the PLD fits. We conclude that the pixel 
mapping technique produces the lowest residual scatter 
for our full-orbit observation data sets, and we therefore 
use this technique in the final version of our analysis. 

In addition to having higher residual scatter than the 


pixel mapping solutions, the PLD fits yield eclipse depth 
and phase curve parameter estimates that often differ 
strongly from the corresponding values derived from fits 
with pixel mapping; these discrepancies sometimes ex¬ 
ceed the 3(7 level. The best-fit values derived from the 
PLD fits also display a higher level of variation across 
different choices of binning, photo metr ic aperture, and 
exponential ramp type (see Section 3.41 than in the case 
of pixel mapping fits. This points toward an inherent in¬ 
stability in the PLD method when fitting full-orbit phase 
curves that may be related to the larger range in star mo¬ 
tions characteristic of such data sets. 

3.4. Exponential ramp correction 

Previous studies using Spitzer/lYlKC have noted a 
short-duration ramp at the beginn ing of each observa¬ 
tion, and again af ter downlinks (e.g. Knutson et al.|20I2 
Lewis et alj20I3 ). The ramp has a characteristic asymp¬ 


totic shape that typically decays to a constant value on 
timescales of an hour or less in the 3.6 and 4.5 /im band- 
passes. We first experimented with removing the first 30 
or 60 minutes of data from each phase curve observation, 
selecting the removal interval that minimizes the residual 
RMS from the best-fit solution binned in five-minute in¬ 
tervals. We find that in both bandpasses, we obtain the 
best results when we do not trim any data from the start 
of the observations. When examining the residual time 
series, we noticed a small ramp visible at the start of the 
3.6 /im observations. We therefore considered whether or 
not our fits might be further improved by the addition 
of an exponential function. 

We experimented with including an exponential ramp 
in our phase curv e model, using the formulation given in 
Agol et al. (2010): 


F = 1 ± aie"‘/“" ± 


Use 


— t/a4 


(7) 


where t is the time since the beginning of the observa¬ 
tion, and ai —04 are correction coefficients. To determine 
whether this function is necessary and if so, how many 
exponential terms to include in the ramp model, we use 
the Bayesian Information Criterion (BIG), defined as 


BlC = x^ + k\nN, 


( 8 ) 


where k is the number of free parameters in the fit, and 
N is the number of data points. By minimizing the BIG, 
we select the type of ramp model that yields the smallest 
residuals without “over-fitting” the data. For the 3.6 /iin 
data, we find that using a single exponential ramp gives 
a marginally lower BIG compared to the no-ramp case, 
while for the 4.5 /tm data, no ramp is needed at all. 
The residuals from the best-fit full phase curve solution, 
shown in Figures H and do not appear to display any 
uncorrected rampdike b3raviors. 

3.5. Parameter fits 

We use a Levenberg-Marquardt least-squares algo¬ 
rithm to fit each full-orbit photometric series to our total 
model light curve, with the intrapixel sensitivity correc¬ 
tion calculated via pixel mapping. In the final version 
of these global fits, we use the updated values for e and 
u! o btain ed from our radial velocity analysis (see Sec¬ 
tion 4.1) as well as the updated orbital period calculated 
from ntting all published transits. The best-fit transit, 
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TABLE 1 

Best-fit Parameters 


Parameter 

3.6 pm 

4.5 pm 

Transit Parameters 

Rp/R* 
tr (BJD)’^ 

0 0941 

U.Uy^iO_Q QQQgg 

2456034.21228t!);)!i()(26 

U.Uy^Zl_Q QQQ7Q 

2456043.18707li(;)ji()!25 

Eclipse Parameters 

1 ®* eclipse depth, di (%) 

tEi (BJD)"’b 

2 "^^ eclipse depth, c ?2 {%)^ 

tE2 (BJD)"’b.t 

u.iooy_Q Q]^Qg 

2456033.05277t)|;(!;('j’))i 

A 1 qqq+0.0060 

U.i?55y_Q QQ4g 

2456035.29948^°;°°°®^ 

0 21 1 c;+0-0185 

U.ZllO_Q Q114 

2456042.02887li(;(ji((!^J 

A nq/^7+0.0096 
U.ZOD/_q Q142 

2456044.274001°;°°°^^ 

Orbital Parameters 

Inclination, i (°) 

Scaled semi-major axis, a/i?» 

84.651°;°® 

6 - 011 !];}^ 

84.611°;°° 

5.98 ±0.13 

Phase Curve Parameters 
Amplitude, ci (xlO“^) 

Phase shift, C 2 (°) 

Maximum flux offset (h)'’ 
Minimum flux offset (h)'’ 

9.701°;°° 

5-711^ 

-1.43 ±0.21 

^•'JO_0.39 

^5D_o.24 

9-4lli 
- 1.01 ± 0.21 

Ramp Parameters 
ai (xlO-4) 

02 (d) 

—4 7+1-8 
‘ -2.0 

U.iUO_Q Q54 

— 


^All tim es are listed in BJDu tc for consistency with other studies; to convert to BJDtt, add 
66.184 s ( Eastman et al.|2010 |. 

*’The center of secondary eclipse times are not corrected for the light travel time across the 
system, At = 35.9 s. 

‘^The maximum and minimum flux offsets are measured relative to the center of secondary eclipse 
time and center of transit time, respectively, and are derived from the phase curve fit parameters 
Cl and C 2 . The maximum flux offset reported is the error-weighted mean of the flux offsets relative 
to the first and second secondary eclipses. 

^These values are computed from fitting the second 3.6 /im secondary eclipse separately using the 
pixel-level decorrelation (PLD) method. This was done in order to remove an anomalous signal 
that occurs during the eclipse in the global 3.6 pm phase curve fit using pixel mapping. 


secondary eclipse, and phase parameters are listed in Ta- 
blefflalong with their uncertainties. Figures]^ andshow 
the Tull-orbit data in the 3.6 and 4.5 pm bandpasses, re¬ 


spectively, with instrumental variations removed. The 
individual eclipse and transit light curves are shown in 
Figures 1^ and 

In the 3.6 pm light curve data, there is an anomalous 
signal in the residuals from the global phase curve fit that 
occurs during the last secondary eclipse. This anomaly 
is characterized by a short ^ 20 -minute dip in the middle 
of the eclipse and was not removed by the pixel map¬ 
ping method, resulting in a S.lcr discrepancy between 
the eclipse depths from the first and second events. Sim¬ 
ilar short-duration anomalies (both positive and nega¬ 
tive) have been reported in Spitzer 3.6 pm data and are 
usually attribu ted to variations in t he width of the stel¬ 
lar target (e.g., Lanotte et al.||2014). Our pixel mapping 


technique accounts for these variations by incorporating 
the noise-pixel parameter pixel sensitivity calculation in 
Eq. (|^, so we conclude that the anomaly in our data is 
likely attributable to some other instrumental effect. 

In order to recover the second 3.6 pm eclipse depth, we 
experimented with fitting the eclipse event separately. 
Selecting a short ^0.2-day segment of the phase curve 
observation surrounding the eclipse, we fit the data to a 
simplified secondary eclipse light curve model using both 
pixel mapping and PLD. This model has three free pa¬ 
rameters — the center of eclipse time tE 2 , the eclipse 
depth d 2 , and a linear slope Ct, which accounts for the 
out-of-eclipse variation in the planet’s brightness. The 
inclination, scaled semi-major axis, and planet-star ra¬ 
dius ratio are fixed at the values derived from the global 
3.6 pm phase curve fit and listed in Tablej^ In these fits, 
we use the same choice of aperture as in the full phase 
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curve fit; in the case of PLD, we optimize for the bin size 
based on the residual scatter and find that a bin size of 
128 yields the lowest residual RMS. Comparing the re¬ 
sults of our fits using pixel mapping and PLD, we see that 
while the residuals from the best-fit solution with pixel 
mapping show a significant anomalous signal similar to 
the one present in the global 3.6 p,m phase curve fit, no 
residual anomaly is evident in the best-fit solution with 
PLD. The range of star positions during this segment of 
data is much less than 0.2 pixels, so we expect PLD to 
perform optimally. The best-fit parameter values from 
the PLD fit are tE 2 = 2456035.(BJDutc), 
d 2 = 0.1894t°;°°59 and c* = 7.7^®;® x lO"-* d-h The 
data with instrumental effects removed and the best-fit 
secondary eclipse solution are shown in Figure ^ We 
note that the best-fit eclipse depth from our individual 
fit of the second 3.6 /rm eclipse is consistent with the 
depth of the first 3.6 ^m eclipse derived from the global 
phase curve fit at the 0.15cr level. Henceforth, we use 
the values from our individual secondary eclipse fit with 
PLD for the center of eclipse and eclipse depth of the 
second 3.6 /im eclipse. 

Taking the error-weighted average of the eclipse depths 
listed in Table P at each wavelength, we arrive at 
0.1882% ± 0.004^o and 0.2247% ± 0.0086% for the 3.6 
and 4.5 /rm bandpasses, respectivel y. These values are 
consistent with the ones reported in Blecic et al. (20131 
in their analysis of previous Spitzer secondary eclipse ob- 
servations to better than la (0.19% ± 0.01% at 3.6 pm 
and 0.224% ± 0.018% at 4.5 ^m). 

We estimate the uncertainties in our best-fit parame¬ 
ters in two wa ys. The first appro ach is the “prayer-bead” 
(PB) method (Gillon et al. 2009), which gives an estimate 
of the contribution of time-correlated noise to the uncer¬ 
tainty. This method entails extracting the residuals from 
the best-fit solution, dividing the residuals into segments, 
and cyclically permuting the residual series segment by 
segment, each time adding the new residual series back 
to the best-fit solution and recomputing the parameters 
using the least-squares algorithm. For each free param¬ 
eter, we create a histogram of the best-fit values from 
every permutation and calculate the uncertainties based 
on the ItT upper and lower bounds from the median. 

The second approach is a Markov chain Monte Carlo 
(MCMC) routine with 10® steps, where we initiate each 
chain at the best-fit solution from the least-squares anal¬ 
ysis. The uncertainty on individual data points is set 
to be the standard deviation of the residuals from the 
best-fit solution. We discard an initial burn-in on each 
chain of length equal to 20% of the chain length, which 
we found ensured the removal of any initial transient be¬ 
havior in a chain, regardless of the choice of initial state. 
The distribution of values for each parameter is close to 
Gaussian, and there are no significant correlations be¬ 
tween pairs of parameters. As in the PB method, we 
set the uncertainties in the fitted parameters to be the 
la upper and lower bounds from the median. For each 
parameter, we choose the larger of the two errors and 
report it in Table We find that the PB errors are con- 
er and r 


sistently larger 


range between 1.0 and 2.7 times that 


of the corresponding MCMC errors. 

The RMS scatter in the best-fit residual series exceeds 
the predicted photon noise limit by a factor of 1.16 at 



-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 

Orbital phase 


Fig. 4.— Top panel: Final 3.6 /im photometric series with instru¬ 
mental variations removed, binned in five-minute intervals (black 
dots). The best-fit total phase, transit, and eclipse light curve is 
overplotted in red. Middle panel: The same data as the upper 
panel, but with an expanded y axis for a clearer view of the phase 
curve. Bottom panel: The residuals from the best-fit solution. 

3.6 pm and 1.14 at 4.5 pm. We estimate the level of red 
noise by calculating the standard deviation of the best-fit 
residuals for various bin sizes, shown in Figure along 
with the inverse square-root dependence of white noise on 
bin size for comparison. On timescales relevant for the 
eclipses and transits (e.g., the ingress/egress timescale 
— ~15 minutes), the red noise increases the RMS by a 
factor of approximately 1.8 at 3.6 ^m and approximately 
1.4 at 4.5 pm. 


4. DISCUSSION 

4.1. Orbital parameters and ephemeris 

We combine the two transit times calculated from our 
phase curve observations with all other published values 


(Joshi et al. 2009 Blecic et al. 2013 Raetz et al.|2015l 


t al.l 

BPd 


to arrive at an updated ephemeris tor the WASP-14 sys- 
tem. Here, we define the zeroth epoch as that of the 
transit nearest in time to the error-weighted mean of all 
measured transit times. The transit observations span 
more than 5.3 years, and by fitting a line through the 
transit times, we derive new, more precise estimates of 
the orbital period P and mid-transit time Tcp: 


fP = 2.24376507 ± 0.00000046 days , . 

= 2455605.65348 ±0.00011 (BJDtdb) ■ 

Figure [T0| shows the observed minus calculated transit 
times derived from these updated ephemeris values. 

We use the secondary eclipse times to obtain a second, 
independent estimate of the orbital period. Carrying out 
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Fig. 5.— Final 4.5 fim photometric series with instrumental vari¬ 
ations removed and the corresponding residual series; see Figure]^ 
for a complete description. 



-0.10 -0.05 0.00 0.05 -0.10 -0.05 0.00 0.05 


Time from Secondary Eclipse Center [d] 


Fig. 6.— Best-fit eclipse light curve data in the 3.6 fim (a—d) 
and 4.5 ^m (e—h) bands after correcting for intrapixel sensitivity 
variations, binned in five-minute intervals (black dots). The best- 
fit model light curves are overplotted in red. The residuals from 
the best-fit solution (c—d, g—h) are shown directly below the cor¬ 
responding light curve data (a—b; e—f). In panels b and d, the 
anomalous residual signal (see text) is highlighted by the gray box. 


a linear fit through the four secondary eclipse times cal¬ 
culated from our phase cu rve data and the two secondary 
eclipse times published in Blecic et al. (20131), we arrive 


at a best-fit value of P = 2.2437060 ± 0.0000017 days. 
This period is consistent with the best-fit transit pe- 



-0.10 -0.05 0.00 0.05 0.10 

Time from Transit Center [d] 


Fig. 7.— Best-fit transit light curve data in the 3.6 /im and 
4.5 fim bands after correcting for intrapixel sensitivity variations, 
binned in five-minute intervals (black dots). The best-fit model 
light curves are overplotted in red. The residuals from the best- 
fit solution are shown directly below the corresponding light curve 
data. 
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Time from Secondary Eclipse Center [d] 


Fig. 8.— Best-fit secondary eclipse light curve data for the sec¬ 
ond 3.6 fim eclipse, derived from fitting the eclipse event separately 
using PLD (black dots). The data has been corrected for intrapixel 
sensitivity variations and is binned in five-minute intervals. The 
best-fit secondary eclipse model light curve is overplotted in red. 
The residuals from the best-fit solution are shown in the bottom 
panel. The gray box indicates the previous location of the anoma¬ 
lous signal (see text), which has been removed by the PLD method. 

riod at the 0.5a level. In Figure [TT] we plot the orbital 
phase of secondary eclipse for aUpublished secondary 
eclipse times using the updated transit ephemeris values 
in Eq. The error-weighted mean orbital phase of 
secondary eclipse is 0.48412 ± 0.00013. 

TABLE 2 

Results from radial velocity fit 

WITH PRIORS ON TRANSIT EPHEMERIS AND ECLIPSE TIMES 


Parameter Value 

Units 

RV Model Parameters 
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Bin size [min] 


Fig. 9.— Plot of the standard deviation of the residuals versus 
bin size for the 3.6 (top panel) and 4.5 /im (bottom panel) data 
sets after removing intrapixel sensitivity effects and dividing out 
the best-fit total light curve solutions (black lines). For comparison, 
the ^1 y/n dependence of white noise on bin size is shown by the red 
lines; the white noise trends are normalized to match the expected 
photon noise limit corresponding to the median photon count over 
each observation. 


TABLE 2 — Continued 




I 

‘■■■1 

- t 

' t 



-400 -300 -200 -100 0 100 200 300 400 

Orbit number 


Fig. 10.— Observed minus calculated transit times for all pub- 
lished observations (red circles are previously-publish ed values; 
|Joshi et al.||2009| [Blecic et al.||2013| | Raetz et al.|[20T^ using the 
updated epnemeris calculated in Section |4.l| me black diamonds 
are the two transit times measured from our phase curve data. 



--i-i 


1 .’ 


-300 -250 -200 -150 190 192 194 196 198 

Orbit number 


Parameter 

Value 

Units 

Pb 

2.24376524 ±4.4E - 07 

days 

Tc,b 

2456034.21261 ±0.00015 

BJDtdb 

Sb 

n 08*^0 +0 0029 
U.USdU _0.0030 


Ub 

252.67 7 - 7 ° 

degrees 

Kb 

986.4 ^2 5 

m s“^ 

7i 

155.8 ±2.6 

m s~^ 

72 

-70.0 tti 

m s~^ 

73 

-73.7 ± 6.8 

m s~^ 

74 

187.4 tl'l 

m s~^ 

7 

0.0025 ±0.0032 

m s~^day~^ 

RV-derived Parameters 


e cosLJ 

— 0 02474 +0-00078 
_o.00074 


esinuj 

—0 0792 +0 0031 
U.U/yz _o.o029 



NOTE — Radial velocity zero point offsets ( 71 — 4 ) derived from 
four separate RV data sets: 1 - Keck/HIRES ijKnutson et al. 


of other variables. 


2014 I, 2 - FIES l |Joshi et al.||2009|l, 3 - SOPHIE < |Joshi et al. 
2009 4 - SOPHIE <|Husnoo et al.|2015|l. See text for description 


Fig. 11.— Orbital phase of secondary eclipse for all publis hed ob¬ 
servations fred circles are previously-published values from [B lecic I 
|et al.|2013|| using the updated ephemeris calculated in Sectionl4.i| 
me black diamonds are the four secondary eclipse times measured 
from our phase curve data. The solid and dashed lines indicate the 
error-weighted mean phase value and corresponding Icr confidence 
bounds, respectively. The horizontal axis has been condensed for 
clarity. 


By combining the updated transit ephemeris and sec¬ 
ondary eclipse times derived from the global phase-curve 
fits with the radial v elocity measurements analyzed in 


Knutson et al. (2014), we can obtain new estimates of 
the orbital eccentricity and pericenter longitude: e = 
O.OSSOto^ooso ^ — 252.67to;77 degre es. These val¬ 
ues ar e consistent with those reported in 
(2014). We use the updated values of or 


Knutson et al.| 
jital eccentric- 


ity andpericenter longitude in our final phase curve fits 
(Table f^. Results from our radial velocity fits are shown 
in Table and Figure [T^ These fits provide new esti- 
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Fig. 12.— Top panel: Phased radial velocity curve of all pub¬ 
lished radial velocity measurements of WASP-14. Bottom panel: 
Corresponding residuals after the radial velocity solution for the 
transiting hot Jupiter is removed. There is no significant linear 
acceleration detected in the data. 


TABLE 3 

Updated Planetary Parameters 


Parameter 

Value 

-Rp j 

0.09419 ± 0.00043 

a/R^, 

5.99 ±0.09 

z n 

84.63 ±0.24 

e 

0.0830;°:°°29 

UJ (°) 

252.67ti(;^? 

Mp (Mjup) 

7.76 ±0.47 

Rp (Rjup) 

1.221 ±0.041 

Pp (g cm-3) 

5.29 ±0.62 

gp (m s- 2 ) 

129 ± 12 

a (AU) 

0.0371 ±0.0011 


mates of the orbital period (Pf,) and center of transit 
time (Tc.b), as well as the semi-amplitude of the planet’s 
radial velocity (Kb), the radial velocity zero point off¬ 
sets for data collected by each of the different spectro¬ 
graphs from which radial velocity measurements of the 
system were obtained ( 71 - 4 ), and the slope ( 7 ) of the 
best-fit radial velocity acceleration. The radial velocity 
slope is consistent with zero, indicating no evidence for 
additional planets in the WASP-14 system. 

From the RV fits we also arrive at updated values of the 
orbital semi-major axis and planet mass: a = 0.0371 ± 
0.0011 AU and Mp = 7.76 ±0.47 Mjup. Using the error- 
weighted best-fit values of a/i?* and Rp/R* from both 
bandpasses, we obtain a new estimate of the planet’s 
radius: Rp = 1.221 ± 0.041 i?jup. A full list of updated 
planetary parameters is give in Table 

4.2. Phase curve fits 

In this section, we combine the results of our global 
fits with model-generated spectra and light curves to 
provide constraints on the atmospheric properties of the 
planet. To compute the relative planet-star flux ratio 


phase curve Fp/F^, we subtract the secondary eclipse 
depth from the best-fit phase curve and divide by the re¬ 
maining flux measured at the center of secondary eclipse, 
which represents the brightness of the star alone. For the 
secondary eclipse depths in each bandpass, we take the 
error-weighted average of the depths listed in Table 
The resulting phase curves are shown in Figure [T^ along 
with the corresponding ±lcr brightness bounds. 

TABLE 4 

Phase Curve Comparison 


Source 

3.6 pva 

4.5 ^m 

Maximum flux ratio [%j 


Measured 

n 1 ©77+0.0094 
U.l5<'/_Q 0;^Qg 

u.zz^y_Q Q090 

Eq. chem. model 

0.2124 

0.2600 

No TiO model 

0.1461 

0.1667 

Minimum flux ratio [%] 


Measured 

< O.OITS’" 

U.UD10_Q QQyg 

Eq. chem. model 

0.0178 

0.0318 

No TiO model 

0.0403 

0.0512 

Phase curve amplitude 1%]^ 


Measured 

> 0.1702=- 

0 1574+00044 

U.1014_Q QQ4g 

Eq. chem. model 

0.1946 

0.2282 

No TiO model 

0.1058 

0.1155 

Maximum flux offset [h]^ 


Measured 

-1.43 ±0.21 

-1.01 ± 0.21 

Eq. chem. model 

±0.3 

±0.2 

No TiO model 

-0.9 

-1.0 

Minimum flux offset [h]^ 


Measured 

—2 

z.uo_Q 30 

— 1 SQ+O'^^ 
4.oy_o 40 

Eq. chem. model 

-3.7 

-2.9 

No TiO model 

-6.2 

-6.2 


^ Values based off of 2cr upper limit of the flux ratio minimum. 
Difference between maximum and minimum flux ratios. 

The maximum and minimum flux offsets are measured relative 


to the center of secondary eclipse time and center of transit time, 
respectively, and are derived from the phase curve fit parameters. 
Negative time offsets of the maximum or minimum flux indicate 
an eastward shift in the location of the hot or cold region in the 
planet’s atmosphere. 


Some of the main quantitative characteristics of the 
phase curves are summarized in Table Notably, the 
fluxes at both Spitzer wavebands peak significantly be¬ 
fore secondary eclipse, implying that the hottest regions 
are shifted eastward from the substellar point. This be¬ 
havior appears to be common on hot Jupiters and, in 


on HD 189733b (IKnutson et al. 

2007 

2009al 2012), HD 

209458b ( 

Zelfem et 

IT pun n 

TKT- 

FT51 

Lewis et al. 

20131), WASP-43b (Stevenson et al. 

2014) 

WASP-12b 

(^Cowan et al. 

2012 

), and Ups And 


Crosshefd et al. 


(Showman & Guillot 

2002 Cooper & Showman 

2005) 

and has now been reproduced in a 

wide varieP 

/ ot gen- 

eral circulation models (GCMs, e.g.. 

Heng et al. 

2011 b|a 

Perna et al.||2012 Lewis et al.||2010t Rauscher ^ 

z IVlenou 
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Fig. 13. — (a) The best-fit 3.6 /im planet-star flux ratio phase curve and Icr brightness bounds (solid and dotted lines). Predicted phase 
curves from the equilibrium chemistry and no TiO SPARC models are plotted with blue circles and red diamonds, respectively, (b) Same 
as (a), but at 4.5 pm. 


2010[ |2012[ 120131 IDobbs-Dixon fc Agol||2013l | Showman] 
et al.|200'^ 2015p. In these models, the eastward hotspot 


displacement results from eastward advection due to a 
fast, eastward-flowing jet stream centered at the equa¬ 
tor. Given such a jet, a significant hot spot offset occurs 
under conditions when the radiative timescale at the pho¬ 
tosphere is comparable to timescales for air to advect hor¬ 
izontally over a planetary radius. The eastward offsets in 
our observations thus provide evidence that WASP-14b 
exhibits an eastward equatorial jet stream. 

It is interesting to compare our observed day-night 
flux differences with those of ot her planets observed to 
date . As mentioned previously , Cow an & Agol (20111 
and Perez-Becker & Showman (2U13) interred that, m 
general, planets that receive higher stellar fluxes exhibit 
larger fractional day-night temperature differences and 
less efficient day-night heat redistribution than planets 
that receive lower stellar fluxes. Averaging the phase 
curve amplitudes at 3.6 and 4.5 ^m (listed in Table®, 
we find that WASP-14b is roughly consistent with tnis 
trend, with a day-night flux difference that is smaller 
than those of highly-irradiated planets like WASP-12b 
and WASP-18b, but larger than those of more weakly- 
irradiated planets like HD 189733b and HD 209458b. 


4.3. Brightness temperature 


We consider four types of atmosphere models. First, 
we use an interpolat ed PHOENIX spec trum for the 


host star WASP-14 (Husser et al. 

2013 

) to calculate 

the brightness temperatures of WASJ 

F-14b in each band 


from the retrieved secondary eclipse depths. The best- 
fit brightness temperatures are 24051*130 K at 3.6 pui 
and 23 931*150 ^ 4.5 /rm. We find that the 3.6 and 
4.5 pm eclipse depths are consistent with a single black- 
body temperature and derive an effective temperature 
of Teff = 2402 ± 35 K from a simultaneous fit to both 
bandpasses. The predicted equilibrium temperature for 
WASP-14b assuming zero albedo is 2220 K if incident 
energy is re-radiated from the dayside only and 1870 K 
if the planet re-radiates the absorbed flux uniformly over 


its entire surface. We can also compare the computed 
brightness temperatures to the effective temperature of 
the dayside in the no-albedo, no-circulation limit assum¬ 
ing each region is a blackbody loc ally in equilibrium with 
the incident stellar flux: 2390 K (|Cowan & Agol|||2011|). 
The high brightness temperatures from the blackbody 
fits therefore suggest that WASP-14b has a very hot day- 
side atmosphere and inefficient day-night recirculation. 


4.4. Dynamical models 

Next, we compare our phase curves and emission 
spectra to theoretical models generated from the three- 
dimensional Substellar and Planetary Atmospheric Radi¬ 
ation and Circulation (SPARC) model to investigate the 
global circulation of WASP-14b. The SPARC model was 
specifically developed with t he study of extrasolar plan¬ 
etary atmospheres in mind (Showm an et al.||2009|). The 
SPARC model couples the MlTgcm ( Adcrott et al. |2UL)4|) 
with the non-g ray radiative transfer model ot |iVlarley &:| 
McKay (19991 in order to self-consistently calculate the 
amount ol heating/cooling at each grid point. In this 
way, the SPARC model does not require advective tun¬ 
ing parameters often employed in one-dimensional radia¬ 
tive transfer models or prescribed pressure-temperature 
profiles utilized in Newtonian cooling schemes employed 
in other circulation models. Our use of an atmospheric 
model that fully couples radiative and dynamical pro¬ 
cesses is especially important for planets on eccentric or¬ 
bits, suc h as WASP-14b, which experience time varia ble 
heating (|Lewis et al.||2010| |2014| |Kataria et al.||2013 1. 

The Sparc models of W ASP-14b presented h ere uti- 
lize the cubed-sphere grid (|Adcroft et al.||2004|) with a 
horizontal grid resolution of roughly 5.625"" in latitude 
and longitude (a so-called low-resolution C16 grid). In 
the vertical direction, the models span pressures ranging 
from 200 bar to 0.2 mbar with 39 layers evenly spaced 
in logp and a top layer that extends to zero pressure. 
Here we consider atmospheric compositions in thermo¬ 
chemical equilibrium both with and without the incor¬ 
poration of the strong optical absorbers TiO and VO 
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Fig. 14.— Comparison of the error-weighted average broadband 
planet-star ratio measured at 3.6, 4.5, and 8.0 iim. (filled black cir¬ 
cles) with SPARC model emission spectra at the time of secondary 
eclipse and the time of transit, corresponding to the dayside and 
nightside of the planet. For the measured nightside 3.6 fim emis¬ 
sion, the 2cr upper limit is shown. Solid lines indicate the predicted 
spectra for the equilibrium chemistry (blue) and no TiO models 
(red). Band-averaged fluxes are overplotted as filled points of the 
same color. The black lines at the bottom represent the photo¬ 
metric band transmission profiles, in arbitrary units. The mea¬ 
sured dayside 3.6 and 4.5 fim planetary fluxes are bounded by the 
two models, suggesting a possible sub-solar abundance of TiO/VO; 
neither model reproduces the low measured 8.0 ^m dayside emis¬ 
sion, which is derived fr om a single seconda ry eclipse measurement 
previously published in|Blecic et al.|(|^01^. Meanwhile, the mea¬ 
sured nightside planetary Iluxes m both the 3.6 and 4.5 /im bands 
are highly discrepant from the model-predicted values, which may 
point toward an enhanced atmospheric C/O ratio. 


(hereafter, “equilibrium chemistry” and “no TiO”, re¬ 
spectively). The presence of TiO/VO in our models 
allows for the develop ment of a dayside inversion layer 
( Fortney et al.|[^08a ). All models adopt solar elemental 
ratios ot heavy elements. 

We assume that WASP-14b is in a pseudo- synchronous 
rotation state (Pr-ot 2.14 days; Hut 1981). As in the 
case of GJ 436b, WASP-14b’s relatively low eccentric¬ 
ity means that assuming pseudo-synchronous rotation 
instead of synchronous rotation is not likely to st rongly 


affect the g lobal circulation patterns that develop (Lewis 

eTalllM^ . 


Fig. 15. — Plot of the dayside and nightside temperature-pressure 
profiles computed by the SPA RC models that correspond to the 
emission spectra in Figure [T^ 


Figure compares the best-fit phase curve derived 
from the Spitzer data in each bandpass with the band- 
averaged light curves generated from the S PARC model 
using the methods of Fortney et al. (20061. Model day- 
side and nightside spectra at the center of secondary 
ecli pse and center of transit times are shown in Fig¬ 
ure along with the corresponding measured flux ratios. 
The corresponding dayside and nightside temperature- 
pressu re p rofiles generated by the models are shown in 
Figure |15| For the dayside planetary emission, we com¬ 
bine the 3.6 and 4.5 pm flux ratios reported i n the present 
work and previously in Blecic et al. (2013) to arrive at 
the error-weighted average values: 0.1^6% ±0.0043% at 
3.6 ^m and 0.2245% ± 0.0078% at 4.5 ^m. We also in¬ 
clude HiejneasuredlftOum flux ratio (0.181% ± 0.021%) 
from Blecic et al. (2013). Table |4| compares the model- 
predicted maximuin and minimum flux ratios and time 
offsets with the values derived from the data. 

In both the 3.6 and 4.5 pm bandpasses, the measured 
dayside planetary emission lies between the equilibrium 
chemistry and no TiO models. A possible explanation is 
that the dayside atmosphere of WASP-14b may contain 
a sub-solar abundance of TiO/VO (possibly due to cold 
trapping), which would result in a weaker temperature 
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inversion than the one predicted b y the solar abundance 
equilibrium chemi stry model (e.g., ^£faigy_gtalj[2008a 


Spiegel et al.|2009] Madhusudhan et al.| 20 ll[ ParmentiCT 


et al. 2013p. Meanwhile, the measured 8.0 ^m dayside 

brightness is not reproduced by either of the models. 

We report offsets in the time of maximum and mini¬ 
mum flux relative to the center of transit and center of 
secondary eclipse times, respectively, as derived from our 
phase curve analysis and the models (Table Negative 
time offsets of the maximum or minimum nux indicate 
an eastward shift in the location of the hot or cold region 
in the planet’s atmosphere. The no TiO model greatly 
overestimates the magnitude of the minimum flux off¬ 
sets while giving good agreement with the maximum flux 
offsets at both wavelengths. The equilibrium chemistry 
model yields a less severe overestimation of the minimum 
flux offset, while predicting a maximum flux offset that is 
opposite the one measured from the data. We find that 
overall the equilibrium chemistry model yields the closer 
match. The no TiO model predicts a significant asym¬ 
metry in the phase curve around the time of minimum 
flux, which is not seen in the Spitzer data. 

The most notable discrepancy between the best-fit 
phase curves and the SPARC model results is in the 
nightside planetary flux. Both the equilibrium chem¬ 
istry and the no TiO model overestimate the planet’s 
nightside brightness at 3.6 ^m, while underestimating 
the brightness at 4.5 fim. The planet’s very low 3.6 fim 
nightside emission indicates a higher atmospheric opac¬ 
ity at that wavelength than is predicted by the models. 
Meanwhile, the higher-than-predicted 4.5 /rm nightside 
emission points toward a slight reduction in the atmo¬ 
spheric opacity at that wavelength relative to the mod¬ 
els. One possible expla nation for both of t hese trends is 
an increased C/0 ratio. Moses et al. (2013) demonstrate 
that increasing the C/0 ratio above equilibrium values 
leads to an excess of CH 4 and a depletion of CO. This 
enhances the opacity in the 3.6 /rm bandpass where CH 4 
has strong vibrational bands, while reducing the opacity 
in the 4.5 ^m bandpass where CO has strong vibrational 
bands. 

An enhanced C/0 ratio can also be invoked to explain 
the higher amplitude of the 3.6 /rm phase curve as com¬ 
pared with the 4.5 /rm phase curve. In this scenario the 
3.6 /im bandpass would probe lower pressure regions that 
have strong day/night contrasts due to a combination of 
short radiative timescales and the possible presence of a 
dayside temperature inversion. A relative depletion of 
CO would shift the 4.5 fim photosphere to higher pres¬ 
sures, where longer radiative timescales facilitate more 
efficient day-night heat transport, resulting in a reduced 
phase curve amplitude. We note, however, that recent 
spectroscopic surveys of M dwarfs suggest that the oc¬ 
currence of high C /O ratios in stellar atmospheres ap¬ 
pears to be low (see Gaidos|2015 and references therein). 
At the same time, rnodels of disk chemistry predict that 
the C/O ratio of the solids and gas in the disk may vary 


as a function of disk r adius (e.g., Oberg et al. 


2011 


Mad- 


husudhan et al.|[2014). In particular, different snowlines 


of oxygen- and carbon-rich ices are expected to result in 
systematic variations in the C/O ratio of gaseous mate¬ 
rial across the protoplanetary disk. Gas giants that de¬ 
rive most of their atmosphere from the gas disk outside of 


the water snowline but inside of the carbon dioxide snow¬ 
line may therefore end up with high atmospheric G/0 
ratios. We also note that Madhusudhan et al. (2011|) 
showed that planets with high C/O ratios should natu¬ 
rally have less TiO and correspondingly weaker tempera¬ 
ture inversion, consistent with our previous discussion of 
WASP -14b’s dayside emission. In addition, [Blecic et al.| 
(2013) show that atmospheric models with high C/O ra¬ 
tios often predict reduced emission at longer wavelengths 
(6 — 8 /im) and may exp lain the low 8.0 /im planetary 
flux ratio (see Figure 14). 

An alternative explanation for the planet’s low 3.6 /im 
nightside emission is the presence of high-altitude silicate 
clouds on the cold nightside in or above the photosphere. 
The formation of an equilibrium silicate cloud is possible 
at the photospheric pressures (~100 mbar) probed on 
the nightside of WA SP-14b by these observations (e.g. 


Visscher et al. 2010). However, it is expected that the 
presence of a thick cloud would suppress both the 3.6 
and 4.5 fiia nightside fluxes from the planet, and the 
resulting planetary emission spectrum would resemble a 
blackbody with an effective temperature corresponding 
to that of the cloud tops. While the low 3.6 pLm night¬ 
side brightness temperature of 1079 K {2a upper limit) 
is consistent with a cold emitting layer, the two night¬ 
side band fluxes together are not consistent with a single 
blackbody spectrum. 

In this scenario, the high 4.5 /rm nightside emis¬ 
sion (corresponding to a brightness temperature of 
1380t™ K) could be explained by introducing an emit¬ 
ting layer of CO in a warmer th ermosphere that is sit - 
uated above the cold cloud tops ( Koskinen et al.p013 ). 
However, it is unclear how a high-altitude temperature 
inversion might arise on the nightside, as there is no inci¬ 
dent stellar irradiation and the efficiency of recirculation 
from the dayside should be relatively weak at low pres¬ 
sures. Therefore, an enhanced C/O ratio in WASP-14b’s 
atmosphere provides a more straightforward explanation 
for the observed differences between WASP-14b’s 3.6 and 
4.5 /im phase curves. Further observations of WASP-14b, 
possibly with Hubble Space Telescope's Wide Field Cam¬ 
era 3, would provide important additional information to 
constrain WASP-14b’s composition and either support or 
refute the high atmospheric C/O ra tio scenario posited 
here for WASP-14b. We note that Lewis et al. (2014) 
also suggested an enhanced C/O ratio to explain the dif- 
ferences seen between the 3.6 and 4.5 /xm phase-curve 
observations of HAT-P-2b. Additional phase-curve ob¬ 
servations of hot Jupiters may reveal similar trends and 
point to a fundamental piece of physics currently miss¬ 
ing from exoplanet atmospheric models and/or planet 
formation theories. The SPARC model utilized in this 
work does not readily accommodate non-solar C/O ra¬ 
tios. The exploration of the effects of different C/O ratios 
on the atmospheric dynamics of hot Jupiters will be the 
topic of future work. 

We also compare the Spitzer dayside and nightside 
emissi on with ID models generated using the meth¬ 
ods of Burrows et al. (2008) with varying recirculation. 
These models assume local thermodynamic equilibrium, 
solar composition, and a plane-parallel atmosphere. A 
heat sink is included at depth (between 0.003 and 0.6 
bars) to redistribute heat from the dayside to the night- 
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side. These models also incorporate a generalized gray 
absorber at low pressures (0 to 0.03 bars), which is 
parametrized with an absorption coefficient, Ate, with 
units of cm^ g“^. This absorber enhances the opacity 
of the planet’s atmosphere at optical wavelengths, rais¬ 
ing the local atmospheric temperature and producing a 
high-altitude temperature inversion. A second dimen¬ 
sionless parameter is used to specify the efficiency of 
energy redistribution, with = 0.5 indicating complete 
redistribution and = 0 signifying redistribution on 
the dayside only. 

In Figure we show the dayside and nightside spec¬ 
tra for various values of Ke and Pn and compare them 
with the measured relative plan etary brightnesses in the 
two Spitzer bands. Figure 17 shows the corresponding 
dayside and nightside temperature-pressure prohles gen¬ 
erated by the model. In these models the presence of the 
gray absorber does not affect day-night recirculation, so 
the nightside spectra depend only on the value of P„. We 
define the measured dayside brightness in each bandpass 
to be the measured secondary eclipse depth, while for the 
the nightside data points, we used the value of the best- 
fit relative planetary phase curve Fp / F^, calculated at the 
mid-transit time in each bandpass. Looking at the day- 
side spectra, we see that the Ke = 0.2 and Pn = 0.1 model 
spectrum comes closest to matching the measured data 
points at 3.6 and 4.5 /im. This suggests that WASP-14b 
has poor day-night recirculation and a moderate thermal 
inversion in the dayside atmosphere. 

On the nightside, as in the SPARC models, none of the 
models reproduce the low measured 3.6 p,m brightness, 
while the measured 4.5 /rm flux ratio is most consistent 
with the Pn = 0.3 model. As discussed earlier in the 
context of the SPARC models, this discrepancy is con¬ 
sistent with the hypothesis that WASP-14b might have 
an enhanced C/0 ratio. We also note that large devi¬ 
ations from solar-like equilibrium chemical composition 
could yield a larger disparity in the pressures probed on 
the dayside and nightside than is accounted for in our 
models. Specifically, one could be probing a deeper than 
expected pressure level on the nightside, where circu¬ 
lation is efficient (high Pn), while probing a higher-up 
(lower pressure) level on the dayside, where the circula¬ 
tion is less efficient (low P„). This scenario is especially 
probable when there is a dayside temperature inversion, 
as our data suggest. 

4.5. Albedo 

Finally, we compare the albedo and recirculation de¬ 
rived from our best-fit phase curves to other hot Jupiters 
for which full-orbit thermal measurements h ave been ob¬ 
tained. Follow ing the methods described in [Schwartz fc| 
Cowan (20151, we use the measured eclipse depths and 
phase curve amplitudes to calculate the error-weighted 
dayside and nightside brightness temperatures, with cor¬ 
rections for the contamination due to reflected light: 
Td = 2312 ± 35 K and r„ = 1299 ± 77 K. Although these 
methods were developed for planets on circular orbits, 
the low eccentricity of WASP-14b places it in a regime 
for which the model is still a reasonable approximation. 
From these temperatures, we derive a Bond albedo of 
Ab = 0 (<0.08 at Icr) and a day-night heat transport 
efficiency of e = 0.23 ± 0.06, where the latter is de¬ 
fined such that e = 0 means no heat recirculation to the 




Fig. 16.— Comparison of the error-weighted average broadband 
planet-star ratio measured at 3.6, 4.5, and 8.0 fim (filled black 
circles) with one-dim ensional atmosphere model spectra following 
I Burrows et al.| (|2008ll at the time of secondary eclipse and the time 
ot transit, corresponding to the dayside and nightside of the planet. 
For the measured nightside 3.6 fim emission, the 2c7 upper limit is 
shown. Solid colored lines indicate the predicted spectra for vari¬ 
ous choices of the parameters ACg Pn , which represent the abun¬ 
dance of a generalized gray absorber in the upper atmosphere and 
the efficiency of energy redistribution, respectively. Correspond¬ 
ing band-averaged points are overplotted in the same color. The 
black lines at the bottom represent the photometric band trans¬ 
mission profiles, in arbitrary units. The measured dayside 3.6 and 
4.5 /im planetary fluxes are most consistent with the model with 
K^e = 0-2 and Pn = 0.1, indicating that WASP-14b has poor day- 
night recirculation and a moderate thermal inversion in the day- 
side atmosphere; none of the models reproduces the low measured 
8.0 /im dayside emission. On the nightside, the measured night¬ 
side planetary fluxes in both the 3.6 and 4.5 /im bands are not 
well-described by any of the models, which may point toward an 
enhanced atmospheric C/O ratio. 


nightside, and e = 1 indicates complete redistribution. 
Figure [l8| shows the location of WASP-14b in albedo- 
recirculation space along with six other exoplanets with 
measured thermal phase curves. WASP-14b’s low day- 
night heat transport efficiency and high irradiation tem¬ 
perature are consistent with the general observed trend 
that high ly-irradiated hot Jupiters have poor heat recir- 


culation (Cowan fc Agol||2011 Perez-Becker & Showman 


20131. 

Both WASP-14b and WASP-18b have much lower es- 
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Fig. 17.— Plot of the dayside and nightside temperature-pressure 
profiles computed by the one-dimension al a tmosphere model used 
to generate the spectra shown in Figure [T^ 


timated albedos than those of other Jovian mass planets 
with thermal phase curve measurements; they are also 
both significantly more massive (7 — 10 vs. ^ 1 Mjup). 
The low estimated albedos of these massive planets, im¬ 
plied by their high thermal emission, may instead be 
indicative of the detection of some amount of thermal 
radiation from the planet’s interior in addition to the 
re-radiated stellar flux. A different cooling and/or mi¬ 
gration history for these massive planets could result in 
an added contribution to the planetary emission from 
the deep atmosphere; the ag e of the WASP-1 4 system 
is relatively young (<1 Gyr; Joshi et al. 2009), so this 
additional internal flux may be due to significant resid¬ 
ual heat of formation. In addition, the higher internal 
fluxes of these higher-mass planets may support stronger 
magnetic fields than on smaller planets, which may lead 
to stronger Ohmic dissipation in their atmosph eres (e.g.. 


Christensen et al. 2009 Batygin et al. 2013), thereby 
providing a source of additional emission on the dayside. 


5. CONCLUSIONS 

In this paper we present the first phase curve obser¬ 
vations of the eccentric hot Jupiter WASP-14b in the 
Spitzer 3.6 and 4.5 ^m bandpasses. We compare two dif¬ 




Fig. 18.— Comparison of composite Icr Bond albedo and recircu¬ 
lation efficiency confidence regions for planets with ther mal phase 
curve observatio ns, as calculated using the methods of |Schwartz| 
I&; Cowan| ||2015ll . The color of the bounding curves indicate irra- 
diation temperature. WASP-14b has intermediate day-night heat 
transport, like other hot Jupiters of similar temperature. How¬ 
ever, it has a very low Bond albedo, like WASP-18b. Both of these 
planets are significantly more massive than the other planets in 
the sample, suggesting that they may be emitting detectable heat 
from the deep atmosphere or interior processes. 


ferent techniques — pixel mapping and pixel-level decor¬ 
relation — for correcting the intrapixel sensitivity effect 
and find that the pixel mapping method yields lower 
residual scatter from the best-fit solution. This is likely 
due to the relatively large range of star positions on the 
pixel throughout the full-orbit observations that may 
make the intrapixel sensitivity effect less amenable to 
modeling through pixel-level decorrelation. We obtain 
best-fit secondary eclipse depths of 0.1882% ± 0.0048% 
and 0.2247% ± 0.0086% at 3.6 and 4.5 /rm, respectively, 
which are consistent with a single blackbody brightness 
temperature of 2402 ± 35 K. These depths are in good 


agreemen t (wit hin Icr) with the ones reported by Ble- 


cic et al. (2013) in their Spitzer secondary eclipse anal¬ 


ysis. Combining the results of our global phase curve 
fits with previous radial velocity measurements, we de¬ 
rive updated values for orbital inclination, orbital eccen¬ 
tricity, longitude of pericenter, orbital semi-major axis, 
planet radius, and planet mass. We also combine our 
measured transit times with previously-published transit 
times to arrive at a more precise estimate of WASP-14b’s 
orbital period: P = 2.24376507 ± 0.00000046 days. 

Comparison of the measured dayside planetary emis¬ 
sion with spectra gener ated from a one-dimen sional ra¬ 
diative transfer model (Burrows et al. 2008) suggests 
relatively inefficient day-night heat recirculation and a 
moderate dayside temperature inversion. The relatively 
high dayside blackbody temperature provides additional 
support for the idea that the day-night circulation is inef¬ 
ficient. The flux maxima precede the secondary eclipses, 
consistent with other hot Jupiter lightcurves and with 
predictions of general circulation models (GCMs), sug¬ 
gesting the possibility of equatorial superrotation on 
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WASP-14b. We also utilize a three-dimensional GCM 
(Showman et al. 2009) to generate theoretical light curves 
tor an atmosphere iii thermochemical equilibrium both 
with and without a dayside temperature inversion. We 
find that the measured amplitude and location of mini¬ 
mum/maximum flux in both bandpasses are more con¬ 
sistent with predictions from the model light curve with 
a dayside thermal inversion. Meanwhile, the measured 
nightside planetary emission at 3.6 and 4.5 /im is not 
adequately described by either the one-dimensional or 
the three-dimensional models. In particular, the very 
low 3.6 /rm nightside planetary flux indicates a signif¬ 
icantly higher atmospheric opacity at that wavelength 
than is predicted by the models, which may point toward 
an enhanced C/0 ratio. In the context of other plan¬ 
ets with full-orbit thermal measurements, we find that 
WASP -14b fits the general trend that highly-irradiated 
hot Jupiters have poor heat recirculation, while the de¬ 
rived Bond albedo is very small (<0.08 at Icr) and the 
planet’s large mass might indicate that WASP-14b is 
emitting residual heat from its formation. 

The question of whether WASP-14b has an enhanced 
C/0 ratio can be further addressed by obtaining mea¬ 
surements of the host star’s C/0 ratio. We note that 
WASP-14 has a near-solar metallicity ([Fe/H] = —0.18± 
0.08; Torres et al.|2012' ), so from a statistical standpoint, 
it is n ot expected to have a high C/0 ratio ( |Teske et al. 


20141. Future work will explore how the assumed 0/0 


ratio in the atmospheres of hot-Jupiters like WASP-14b 
affects global circulation patterns and day/night temper¬ 
ature contrasts. The discrepancy between the model- 
predicted and measured nightside planetary flux under¬ 
lines the need for further exploration of the available 
parameter space in both ID and GCM models of hot 
Jupiter atmospheres. Assessing the effects of non-solar 


chemistry, and specifically different C/0 ratios, will give 
us new insight into the interplay between various atmo¬ 
spheric properties and the resultant planetary emission. 
These studies promise to greatly enhance the ability of 
numerical models to explain features in the growing body 
of exoplanet phase curves. 

Obtaining more infrared phase curves of massive plan¬ 
ets and calculating their Bond albedo will allow us to 
determine whether low Bond albedo is strongly corre¬ 
lated with high planet mass, as is the case for WASP-14b 
and WASP-18b, and further consider the contribution of 
residual heat from formation in the overall emission of 
massive hot Jupiters. A better understanding of the role 
of residual heat from formation in the flux of massive 
planets is also important for directly-imaged planets, as 
cooling curves are used to estimate the planetary mass 
based on th e age of the system an d the measured lumi¬ 
nosity (e.g., Fortney et al. 2008b). If we can determine 
whether current cooling rnodels can reproduce the ob¬ 
served emission from WASP-14b, it would serve as an in¬ 
dependent confirmation of the evolutionary models used 
in studies of directly imaged planets. 
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